In [ ]:
import os
import matplotlib.pyplot as plt
%matplotlib inline
In [ ]:
from rmgpy.chemkin import *
from rmgpy.species import Species
In [ ]:
from afm.simulator import OdeSimulator
import afm.utils
import afm.simulator
In [ ]:
temperature = 673.15 # K
pressure = 350*3 # bar
initial_mol_fraction = {
"ArCCCCR":1.0,
"LCCCCR":1.0,
"LCCCC":1.0
}
hr = 14
termination_time = 3600*hr # hrs
In [ ]:
model = 'two-sided_newcut1'
working_dir = os.path.join('../', 'data', 'pdd_chemistry', model)
chemkin_path = os.path.join(working_dir, 'chem_annotated.inp')
species_dict_path = os.path.join(working_dir, 'species_dictionary.txt')
smiles_dict_path = os.path.join(working_dir, 'fragment_smiles.txt')
In [ ]:
ode_simulator = OdeSimulator(chemkin_path,
species_dict_path,
smiles_dict_path,
temperature,
pressure)
In [ ]:
alldata = ode_simulator.simulate(initial_mol_fraction, termination_time)
In [ ]:
results_path = os.path.join(working_dir, 'results')
if not os.path.exists(results_path):
os.mkdir(results_path)
In [ ]:
# prepare moles data
time, dataList, _ = alldata[0]
TData = dataList[0]
PData = dataList[1]
VData = dataList[2]
total_moles = PData.data*VData.data/8.314/TData.data
moles_dict = {}
for data in dataList[3:]:
spe_label = data.label
moles_dict[spe_label] = max(data.data[-1]*total_moles[-1],0)
In [ ]:
ArCCCCR_mf = dataList[3].data
print dataList[3].label
ArCCCCR_moles = ArCCCCR_mf*total_moles
ArCCCCR_conv = (ArCCCCR_moles[0]-ArCCCCR_moles)/ArCCCCR_moles[0]
In [ ]:
plt.plot(time.data, ArCCCCR_conv)
numpy.savetxt(os.path.join(results_path, 'reactant_conv.csv'), (time.data, ArCCCCR_conv))
In [ ]:
from afm.simulator import categorize_fragments
In [ ]:
fragmental_weight_distri = ode_simulator.get_molecular_weight_distribution(alldata)
In [ ]:
mws = [tup[0]*1000 for tup in fragmental_weight_distri]
moles = [tup[1] for tup in fragmental_weight_distri]
molefracs = moles/sum(moles)
In [ ]:
numpy.savetxt(os.path.join(results_path, 'mwd_{0}hr.csv'.format(hr)), (mws, molefracs))
In [ ]: